setwd(dirname(rstudioapi::getActiveDocumentContext()$path))

## Import packages
require(cregg)
require(ggplot2)
require(haven)
require(psych)

#################
## Import Data ##
#################

## Experiment 1
dcj_exp1 <- readRDS("elecvspref_jan22_conjoint_v4.1.rds")

## Experiment 2
dcj_exp2 <- readRDS("elecvspref_mar16_conjoint_v4.1.rds")

#######################################
## Adding Gender Roles View Variable ##
#######################################

## First Experiment Data
d_exp1 <- subset(dcj_exp1, task==1 & profile==1)
### Gender Roles View Variables
round(prop.table(table(d_exp1$genderval_1,useNA="always")),3) 
round(prop.table(table(d_exp1$genderval_2,useNA="always")),3) 
round(prop.table(table(d_exp1$genderval_3,useNA="always")),3) 
round(prop.table(table(d_exp1$genderval_4,useNA="always")),3) 
round(prop.table(table(d_exp1$genderval_5,useNA="always")),3) 
round(prop.table(table(d_exp1$genderval_6,useNA="always")),3) 
round(prop.table(table(d_exp1$genderval_7,useNA="always")),3) 
round(prop.table(table(d_exp1$genderval_8,useNA="always")),3) 
round(prop.table(table(d_exp1$genderval_9,useNA="always")),3) 
genderval_1a <- ifelse(is.na(d_exp1$genderval_1),2.5,d_exp1$genderval_1)
genderval_2a <- ifelse(is.na(d_exp1$genderval_2),2.5,d_exp1$genderval_2)
genderval_3a <- ifelse(is.na(d_exp1$genderval_3),2.5,d_exp1$genderval_3)
genderval_4a <- ifelse(is.na(d_exp1$genderval_4),2.5,d_exp1$genderval_4)
genderval_5a <- ifelse(is.na(d_exp1$genderval_5),2.5,d_exp1$genderval_5)
genderval_6a <- ifelse(is.na(d_exp1$genderval_6),2.5,d_exp1$genderval_6)
genderval_7a <- ifelse(is.na(d_exp1$genderval_7),2.5,d_exp1$genderval_7)
genderval_8a <- ifelse(is.na(d_exp1$genderval_8),2.5,d_exp1$genderval_8)
genderval_9a <- ifelse(is.na(d_exp1$genderval_9),2.5,d_exp1$genderval_9)
### Create the measurement through factor analysis
d_exp1$genderval <- (((5-genderval_1a) + genderval_2a + genderval_3a +
                   (5-genderval_4a) + genderval_5a + genderval_6a +
                   genderval_7a + (5-genderval_8a) + genderval_9a)/9 - 1)/3

## Second Experiment Data
d_exp2 <- subset(dcj_exp2, task==1 & profile==1)
### Gender Roles View Variables
round(prop.table(table(d_exp2$genderval_1,useNA="always")),3) 
round(prop.table(table(d_exp2$genderval_2,useNA="always")),3) 
round(prop.table(table(d_exp2$genderval_3,useNA="always")),3) 
round(prop.table(table(d_exp2$genderval_4,useNA="always")),3) 
round(prop.table(table(d_exp2$genderval_5,useNA="always")),3) 
round(prop.table(table(d_exp2$genderval_6,useNA="always")),3) 
round(prop.table(table(d_exp2$genderval_7,useNA="always")),3) 
round(prop.table(table(d_exp2$genderval_8,useNA="always")),3) 
round(prop.table(table(d_exp2$genderval_9,useNA="always")),3) 
genderval_1b <- ifelse(is.na(d_exp2$genderval_1),2.5,d_exp2$genderval_1)
genderval_2b <- ifelse(is.na(d_exp2$genderval_2),2.5,d_exp2$genderval_2)
genderval_3b <- ifelse(is.na(d_exp2$genderval_3),2.5,d_exp2$genderval_3)
genderval_4b <- ifelse(is.na(d_exp2$genderval_4),2.5,d_exp2$genderval_4)
genderval_5b <- ifelse(is.na(d_exp2$genderval_5),2.5,d_exp2$genderval_5)
genderval_6b <- ifelse(is.na(d_exp2$genderval_6),2.5,d_exp2$genderval_6)
genderval_7b <- ifelse(is.na(d_exp2$genderval_7),2.5,d_exp2$genderval_7)
genderval_8b <- ifelse(is.na(d_exp2$genderval_8),2.5,d_exp2$genderval_8)
genderval_9b <- ifelse(is.na(d_exp2$genderval_9),2.5,d_exp2$genderval_9)
### Create the measurement through factor analysis
d_exp2$genderval <- (((5-genderval_1b) + genderval_2b + genderval_3b +
                        (5-genderval_4b) + genderval_5b + genderval_6b +
                        genderval_7b + (5-genderval_8b) + genderval_9b)/9 - 1)/3

## Factor Analysis
tmp <- fa(cbind(c(genderval_1a,genderval_1b),
                c(genderval_2a,genderval_2b),
                c(genderval_3a,genderval_3b),
                c(genderval_4a,genderval_4b),
                c(genderval_5a,genderval_5b),
                c(genderval_6a,genderval_6b),
                c(genderval_7a,genderval_7b),
                c(genderval_8a,genderval_8b),
                c(genderval_9a,genderval_9b)), 
          scores="Bartlett",fm="pa")
tmp$loadings # as it is

## Check Correlation and Insert Factors Scores
### Experiment 1
plot(d_exp1$genderval, tmp$scores[1:nrow(d_exp1),1])
cor(d_exp1$genderval, tmp$scores[1:nrow(d_exp1),1])
d_exp1$genderval <- as.numeric(tmp$scores[1:nrow(d_exp1),1])
# hist(d_exp1$genderval)
### Experiment 2
plot(d_exp2$genderval, tmp$scores[(nrow(d_exp1)+1):(nrow(d_exp1)+nrow(d_exp2))])
cor(d_exp2$genderval, tmp$scores[(nrow(d_exp1)+1):(nrow(d_exp1)+nrow(d_exp2))])
d_exp2$genderval <- as.numeric(tmp$scores[(nrow(d_exp1)+1):(nrow(d_exp1)+nrow(d_exp2)),1])
# hist(d_exp2$genderval)

## Median Split Variables
#### Splitting Point (pooled measure)
splitpoint <- median(c(d_exp1$genderval,d_exp2$genderval))
### Experiment 1
d_exp1$genderval_2cat <- 
  ifelse(d_exp1$genderval>=splitpoint,2,1)
d_exp1$genderval_2cat <- 
  factor(d_exp1$genderval_2cat, labels = c("Liberal","Traditional"))
round(prop.table(table(d_exp1$genderval_2cat)),3)
### Experiment 2
d_exp2$genderval_2cat <- 
  ifelse(d_exp2$genderval>=splitpoint,2,1)
d_exp2$genderval_2cat <- 
  factor(d_exp2$genderval_2cat, labels = c("Liberal","Traditional"))
round(prop.table(table(d_exp2$genderval_2cat)),3)

## Only keep relevant variables in datasets
### Experiment 1
tmp <- d_exp1[,c("Response.ID","genderval","genderval_2cat")]
rm(d_exp1)
d_exp1 <- tmp
### Experiment 2
tmp <- d_exp2[,c("Response.ID","genderval","genderval_2cat")]
rm(d_exp2)
d_exp2 <- tmp

##################
## Experiment 1 ##
##################

## Clean Data
dim(dcj_exp1)
## Experimental Assignment Mechanically Missing
dcj_exp1 <- subset(dcj_exp1, !is.na(dcj_exp1$candgender))
dim(dcj_exp1)
## Non-Responses
dcj_exp1 <- subset(dcj_exp1, !is.na(dcj_exp1$selected))
dim(dcj_exp1)

## Add Gender Roles View variable
dcj_exp1$genderval <- 
  d_exp1$genderval[match(dcj_exp1$Response.ID,d_exp1$Response.ID)]
dcj_exp1$genderval_2cat <- 
  d_exp1$genderval_2cat[match(dcj_exp1$Response.ID,d_exp1$Response.ID)]

## Adjust Variables
dcj_exp1$dv <- 
  factor(ifelse(dcj_exp1$electability==1, "Expectation", "Preference"),
         levels = c("Preference", "Expectation"))
dcj_exp1$candgender <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candgender),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candgender)))
dcj_exp1$candparty <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candparty),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candparty)))
dcj_exp1$candage <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candage),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candage)))
dcj_exp1$candexperience <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candexperience),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candexperience)))
dcj_exp1$candedu <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candedu),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candedu)))
dcj_exp1$candmarriage <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candmarriage),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candmarriage)))
dcj_exp1$candkidage <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candkidage),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candkidage)))
dcj_exp1$candlivetype <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp1$candlivetype),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp1$candlivetype)))

## Variable Labels 
vlab_exp1 <- list("candgender"="Gender",
             "candparty"="Party Affliation",
             "candage"="Age",
             "candexperience"="Political Experience",
             "candedu"="Educational Attainment",
             "candmarriage"="Marital Status",
             "candkidage"="Youngest Child",
             "candlivetype"="Residence")

## Marginal Means ##

resexp1_each <- cj(selected ~ candgender + candparty + candage + 
                 candexperience + candedu + candmarriage + 
                 candkidage + candlivetype, 
               data = dcj_exp1,
               feature_labels = vlab_exp1,
               by = ~dv+genderval_2cat, 
               id = ~Response.ID, estimate="mm", h0=0.5)
head(resexp1_each)
resexp1_diff1 <- cj(selected ~ candgender + candparty + candage + 
                 candexperience + candedu + candmarriage + 
                 candkidage + candlivetype, 
               data = subset(dcj_exp1, genderval_2cat=="Liberal"),
               feature_labels = vlab_exp1,
               by = ~dv, 
               id = ~Response.ID, estimate="mm_diff")
resexp1_diff1$genderval_2cat <- "Liberal"
head(resexp1_diff1)
resexp1_diff2 <- cj(selected ~ candgender + candparty + candage + 
                      candexperience + candedu + candmarriage + 
                      candkidage + candlivetype, 
                    data = subset(dcj_exp1, genderval_2cat=="Traditional"),
                    feature_labels = vlab_exp1,
                    by = ~dv, 
                    id = ~Response.ID, estimate="mm_diff")
resexp1_diff2$genderval_2cat <- "Traditional"
head(resexp1_diff2)

resexp1_both <- rbind(resexp1_each,resexp1_diff1,resexp1_diff2)
resexp1_both$BY <- 
  factor(gsub("\\*\\*\\*.*$","",resexp1_both$BY), 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp1_both$genderval_2cat <- 
  factor(resexp1_both$genderval_2cat, 
         levels = c("Liberal","Traditional"))
resexp1_both$level <- factor(resexp1_both$level, 
                         levels = rev(levels(resexp1_both$level)))

setfacetexp1_vline <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             VL = rep(c(0.5,0.5,0),each=2))
setfacetexp1_lims <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             level = factor("Male", levels(resexp1_both$level)), 
             feature = factor("Gender", levels(resexp1_both$feature)),
             lims = c(0.3,0.7,0.3,0.7,-0.15,0.15))

## Plot!
## Figure A11 (Appendix D.3)
require(ggplot2)
p_resexp1 <- 
  ggplot(resexp1_both) + 
  geom_vline(data=setfacetexp1_vline, aes(xintercept=VL), 
             linetype=2, color="gray70") + 
  geom_errorbarh(aes(xmin=lower,xmax=upper,y=level, 
                     color=genderval_2cat), height=0, linewidth = 0.5,
                 position = position_dodge(width=-0.8)) +
  geom_point(aes(x=estimate,y=level,
                 shape=genderval_2cat), size=1,
             position = position_dodge(width=-0.8)) + 
  geom_point(data=setfacetexp1_lims, aes(x=lims,y=level), 
             size=1, alpha=0) + 
  scale_shape_discrete(name="Views on Gender Roles") + 
  scale_color_discrete(name="Views on Gender Roles") + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="Marginal Means Estimates (with 95% Confidence Intervals)",
       y=NULL) + 
  theme_bw(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom")
p_resexp1
ggsave("resexp1_genderval_v4.pdf", width=8, height=7)

### Extracting p-values
### Figure A12 (Appendix D.3)
require(ggplot2)
p_resexp1_pval <- 
  ggplot(resexp1_both) + 
  geom_label(aes(x=ifelse(genderval_2cat=="Liberal",0.49,0.51),
                 y=level,label=sprintf("%.03f",p), 
                 fill=genderval_2cat), size=3) + 
  geom_point(aes(x=0.48,y=level), size=1, alpha=0) + 
  geom_point(aes(x=0.52,y=level), size=1, alpha=0) + 
  scale_fill_discrete(name="Views on Gender Roles") + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="p-values from significance test",
       y=NULL,
       caption="Note: H0=0.5 for preference and expectation tasks. H0=0 for preference-expecation gap.") + 
  theme_minimal(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        axis.text.x = element_blank())
p_resexp1_pval
ggsave("resexp1_genderval_pval_v4.pdf", width=8, height=7)

## AMCE (only estimation) ##

resexp1_amce_each <- cj(selected ~ candgender + candparty + candage + 
                     candexperience + candedu + candmarriage + 
                     candkidage + candlivetype, 
                   data = dcj_exp1,
                   feature_labels = vlab_exp1,
                   by = ~dv+genderval_2cat, 
                   id = ~Response.ID, estimate="amce")
head(resexp1_amce_each)
resexp1_amce_diff1 <- cj(selected ~ candgender + candparty + candage + 
                     candexperience + candedu + candmarriage + 
                     candkidage + candlivetype, 
                   data = subset(dcj_exp1, genderval_2cat=="Liberal"),
                   feature_labels = vlab_exp1,
                   by = ~dv, 
                   id = ~Response.ID, estimate="amce_diff")
resexp1_amce_diff1$genderval_2cat <- "Liberal"
head(resexp1_amce_diff1)
resexp1_amce_diff2 <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candmarriage + 
                           candkidage + candlivetype, 
                         data = subset(dcj_exp1, genderval_2cat=="Traditional"),
                         feature_labels = vlab_exp1,
                         by = ~dv, 
                         id = ~Response.ID, estimate="amce_diff")
resexp1_amce_diff2$genderval_2cat <- "Traditional"
head(resexp1_amce_diff2)

resexp1_amce_both <- rbind(resexp1_amce_each,resexp1_amce_diff1,resexp1_amce_diff2)
resexp1_amce_both$BY <- 
  factor(gsub("\\*\\*\\*.*$","",resexp1_amce_both$BY), 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp1_amce_both$genderval_2cat <- 
  factor(resexp1_amce_both$genderval_2cat, 
         levels = c("Liberal","Traditional"))
resexp1_amce_both$level <- factor(resexp1_amce_both$level, 
                             levels = rev(levels(resexp1_amce_both$level)))

##################
## Experiment 2 ##
##################

## Clean Data
dim(dcj_exp2)
## Experimental Assignment Mechanically Missing
dcj_exp2 <- subset(dcj_exp2, !is.na(dcj_exp2$candgender))
dim(dcj_exp2)
## Nonresponses
dcj_exp2 <- subset(dcj_exp2, !is.na(dcj_exp2$selected))
dim(dcj_exp2)

## Add Gender Roles View variable
dcj_exp2$genderval <- 
  d_exp2$genderval[match(dcj_exp2$Response.ID,d_exp2$Response.ID)]
dcj_exp2$genderval_2cat <- 
  d_exp2$genderval_2cat[match(dcj_exp2$Response.ID,d_exp2$Response.ID)]

## Adjust Variables
dcj_exp2$dv <- 
  factor(ifelse(dcj_exp2$electability==1, "Expectation", "Preference"),
         levels = c("Preference", "Expectation"))
dcj_exp2$candgender <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candgender),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candgender)))
dcj_exp2$candparty <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candparty),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candparty)))
dcj_exp2$candage <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candage),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candage)))
dcj_exp2$candexperience <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candexperience),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candexperience)))
dcj_exp2$candedu <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candedu),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candedu)))
dcj_exp2$candfam <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candfam),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candfam)))
dcj_exp2$candlivetype <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candlivetype),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candlivetype)))
dcj_exp2$candorigin <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candorigin),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candorigin)))
dcj_exp2$candimppol <- 
  factor(gsub("^\\(.*\\) ","", dcj_exp2$candimppol),
         levels = gsub("^\\(.*\\) ","", levels(dcj_exp2$candimppol)))

## Variable Labels 
vlab_exp2 <- list("candgender"="Gender",
                  "candparty"="Party Affliation",
                  "candage"="Age",
                  "candexperience"="Political Experience",
                  "candedu"="Educational Attainment",
                  "candfam"="Family",
                  "candorigin"="Native Place",
                  "candlivetype"="Residence",
                  "candimppol"="Policy Focus")

## house (Marginal Means) ##

resexp2_each_house <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="house"),
                         feature_labels = vlab_exp2,
                         by = ~dv+genderval_2cat, 
                         id = ~Response.ID, estimate="mm", h0=0.5)
head(resexp2_each_house)
resexp2_diff_house1 <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="house" & 
                                         genderval_2cat == "Liberal"),
                         feature_labels = vlab_exp2,
                         by = ~dv, 
                         id = ~Response.ID, estimate="mm_diff")
resexp2_diff_house1$genderval_2cat <- "Liberal"
head(resexp2_diff_house1)
resexp2_diff_house2 <- cj(selected ~ candgender + candparty + candage + 
                            candexperience + candedu + candfam + 
                            candlivetype + candorigin + candimppol, 
                          data = subset(dcj_exp2, conjoint_eleclevel=="house" & 
                                          genderval_2cat == "Traditional"),
                          feature_labels = vlab_exp2,
                          by = ~dv, 
                          id = ~Response.ID, estimate="mm_diff")
resexp2_diff_house2$genderval_2cat <- "Traditional"
head(resexp2_diff_house2)

resexp2_both_house <- rbind(resexp2_each_house,
                            resexp2_diff_house1,
                            resexp2_diff_house2)
resexp2_both_house$BY <- 
  factor(gsub("\\*\\*\\*.*$", "", resexp2_both_house$BY), 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp2_both_house$genderval_2cat <- 
  factor(resexp2_both_house$genderval_2cat, 
         levels = c("Liberal","Traditional"))
resexp2_both_house$level <- 
  factor(resexp2_both_house$level, 
         levels = rev(levels(resexp2_both_house$level)))

setfacetexp2_vline_house <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             VL = rep(c(0.5,0.5,0),each=2))
setfacetexp2_lims_house <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             level = factor("Male", levels(resexp2_both_house$level)), 
             feature = factor("Gender", levels(resexp2_both_house$feature)),
             lims = c(0.3,0.7,0.3,0.7,-0.15,0.15))

## Plot!
## Figure A13 (Appendix D.3)
require(ggplot2)
p_resexp2_house <- 
  ggplot(resexp2_both_house) + 
  geom_vline(data=setfacetexp2_vline_house, aes(xintercept=VL), 
             linetype=2, color="gray70") + 
  geom_errorbarh(aes(xmin=lower,xmax=upper,y=level, 
                     color=genderval_2cat), height=0, linewidth = 0.5,
                 position = position_dodge(width=-0.8)) +
  geom_point(aes(x=estimate,y=level,
                 shape=genderval_2cat), size=1,
             position = position_dodge(width=-0.8)) + 
  geom_point(data=setfacetexp2_lims_house, aes(x=lims,y=level), 
             size=1, alpha=0) + 
  scale_shape_discrete(name="Views on Gender Roles") + 
  scale_color_discrete(name="Views on Gender Roles") + 
  # scale_color_manual(name="Women's Issue Minister", values=c("gray50","black",
  #                                                            "orange1","darkorange4")) + 
  # scale_shape_manual(name="Women's Issue Minister", values=c(16,16,17,17)) + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="Marginal Means Estimates (with 95% Confidence Intervals)",
       y=NULL) + 
  theme_bw(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom")
p_resexp2_house
ggsave("resexp2_house_genderval_v4.pdf", width=8, height=7)

### Extracting p-values
### Figure A14 (Appendix D.3)
require(ggplot2)
p_resexp2_house_pval <- 
  ggplot(resexp2_both_house) + 
  geom_label(aes(x=ifelse(genderval_2cat=="Liberal",0.49,0.51),
                y=level,label=sprintf("%.03f",p), 
                fill=genderval_2cat), size=3) + 
  geom_point(aes(x=0.48,y=level), size=1, alpha=0) + 
  geom_point(aes(x=0.52,y=level), size=1, alpha=0) + 
  scale_fill_discrete(name="Views on Gender Roles") + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="p-values from significance test",
       y=NULL,
       caption="Note: H0=0.5 for preference and expectation tasks. H0=0 for preference-expecation gap.") + 
  theme_minimal(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        axis.text.x = element_blank())
p_resexp2_house_pval
ggsave("resexp2_house_genderval_pval_v4.pdf", width=8, height=7.5)

## house (AMCE, only estimation) ##

resexp2_amce_each_house <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="house"),
                         feature_labels = vlab_exp2,
                         by = ~dv+genderval_2cat, 
                         id = ~Response.ID, estimate="amce")
head(resexp2_amce_each_house)
resexp2_amce_diff_house1 <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="house" & 
                                         genderval_2cat=="Liberal"),
                         feature_labels = vlab_exp2,
                         by = ~dv, 
                         id = ~Response.ID, estimate="amce_diff")
resexp2_amce_diff_house1$genderval_2cat <- "Liberal"
head(resexp2_amce_diff_house1)
resexp2_amce_diff_house2 <- cj(selected ~ candgender + candparty + candage + 
                                candexperience + candedu + candfam + 
                                candlivetype + candorigin + candimppol, 
                              data = subset(dcj_exp2, conjoint_eleclevel=="house" & 
                                              genderval_2cat=="Traditional"),
                              feature_labels = vlab_exp2,
                              by = ~dv, 
                              id = ~Response.ID, estimate="amce_diff")
resexp2_amce_diff_house2$genderval_2cat <- "Traditional"
head(resexp2_amce_diff_house2)

resexp2_amce_both_house <- rbind(resexp2_amce_each_house,
                            resexp2_amce_diff_house1,
                            resexp2_amce_diff_house2)
resexp2_amce_both_house$BY <- 
  factor(gsub("\\*\\*\\*.*$","",resexp2_amce_both_house$BY), 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp2_amce_both_house$genderval_2cat <- 
  factor(resexp2_amce_both_house$genderval_2cat, 
         levels = c("Liberal","Traditional"))
resexp2_amce_both_house$level <- 
  factor(resexp2_amce_both_house$level, 
         levels = rev(levels(resexp2_amce_both_house$level)))

## municipal (Marginal Means) ##

resexp2_each_municipal <- cj(selected ~ candgender + candparty + candage + 
                           candexperience + candedu + candfam + 
                           candlivetype + candorigin + candimppol, 
                         data = subset(dcj_exp2, conjoint_eleclevel=="municipal"),
                         feature_labels = vlab_exp2,
                         by = ~dv+genderval_2cat, 
                         id = ~Response.ID, estimate="mm", h0=0.5)
head(resexp2_each_municipal)
resexp2_diff_municipal1 <- cj(selected ~ candgender + candparty + candage + 
                            candexperience + candedu + candfam + 
                            candlivetype + candorigin + candimppol, 
                          data = subset(dcj_exp2, conjoint_eleclevel=="municipal" & 
                                          genderval_2cat == "Liberal"),
                          feature_labels = vlab_exp2,
                          by = ~dv, 
                          id = ~Response.ID, estimate="mm_diff")
resexp2_diff_municipal1$genderval_2cat <- "Liberal"
head(resexp2_diff_municipal1)
resexp2_diff_municipal2 <- cj(selected ~ candgender + candparty + candage + 
                            candexperience + candedu + candfam + 
                            candlivetype + candorigin + candimppol, 
                          data = subset(dcj_exp2, conjoint_eleclevel=="municipal" & 
                                          genderval_2cat == "Traditional"),
                          feature_labels = vlab_exp2,
                          by = ~dv, 
                          id = ~Response.ID, estimate="mm_diff")
resexp2_diff_municipal2$genderval_2cat <- "Traditional"
head(resexp2_diff_municipal2)

resexp2_both_municipal <- rbind(resexp2_each_municipal,
                            resexp2_diff_municipal1,
                            resexp2_diff_municipal2)
resexp2_both_municipal$BY <- 
  factor(gsub("\\*\\*\\*.*$", "", resexp2_both_municipal$BY), 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp2_both_municipal$genderval_2cat <- 
  factor(resexp2_both_municipal$genderval_2cat, 
         levels = c("Liberal","Traditional"))
resexp2_both_municipal$level <- 
  factor(resexp2_both_municipal$level, 
         levels = rev(levels(resexp2_both_municipal$level)))

setfacetexp2_vline_municipal <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             VL = rep(c(0.5,0.5,0),each=2))
setfacetexp2_lims_municipal <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             level = factor("Male", levels(resexp2_both_municipal$level)), 
             feature = factor("Gender", levels(resexp2_both_municipal$feature)),
             lims = c(0.3,0.7,0.3,0.7,-0.15,0.15))

## Plot!
## Figure A15 (Appendix D.3)
require(ggplot2)
p_resexp2_municipal <- 
  ggplot(resexp2_both_municipal) + 
  geom_vline(data=setfacetexp2_vline_municipal, aes(xintercept=VL), 
             linetype=2, color="gray70") + 
  geom_errorbarh(aes(xmin=lower,xmax=upper,y=level, 
                     color=genderval_2cat), height=0, linewidth = 0.5,
                 position = position_dodge(width=-0.8)) +
  geom_point(aes(x=estimate,y=level,
                 shape=genderval_2cat), size=1,
             position = position_dodge(width=-0.8)) + 
  geom_point(data=setfacetexp2_lims_municipal, aes(x=lims,y=level), 
             size=1, alpha=0) + 
  scale_shape_discrete(name="Views on Gender Roles") + 
  scale_color_discrete(name="Views on Gender Roles") + 
  # scale_color_manual(name="Women's Issue Minister", values=c("gray50","black",
  #                                                            "orange1","darkorange4")) + 
  # scale_shape_manual(name="Women's Issue Minister", values=c(16,16,17,17)) + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="Marginal Means Estimates (with 95% Confidence Intervals)",
       y=NULL) + 
  theme_bw(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom")
p_resexp2_municipal
ggsave("resexp2_municipal_genderval_v4.pdf", width=8, height=7)

### Extracting p-values
### Figure A16 (Appendix D.3)
require(ggplot2)
p_resexp2_municipal_pval <- 
  ggplot(resexp2_both_municipal) + 
  geom_label(aes(x=ifelse(genderval_2cat=="Liberal",0.49,0.51),
                 y=level,label=sprintf("%.03f",p), 
                 fill=genderval_2cat), size=3) + 
  geom_point(aes(x=0.48,y=level), size=1, alpha=0) + 
  geom_point(aes(x=0.52,y=level), size=1, alpha=0) + 
  scale_fill_discrete(name="Views on Gender Roles") + 
  facet_grid(feature~BY, scales = "free", space = "free") + 
  labs(x="p-values from significance test",
       y=NULL,
       caption="Note: H0=0.5 for preference and expectation tasks. H0=0 for preference-expecation gap.") + 
  theme_minimal(base_family = "serif") + 
  theme(strip.text.y = element_text(angle=0),
        plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        axis.text.x = element_blank())
p_resexp2_municipal_pval
ggsave("resexp2_municipal_genderval_pval_v4.pdf", width=8, height=7.5)

## municipal (AMCE, only estimation) ##

resexp2_amce_each_municipal <- cj(selected ~ candgender + candparty + candage + 
                                candexperience + candedu + candfam + 
                                candlivetype + candorigin + candimppol, 
                              data = subset(dcj_exp2, conjoint_eleclevel=="municipal"),
                              feature_labels = vlab_exp2,
                              by = ~dv+genderval_2cat, 
                              id = ~Response.ID, estimate="amce")
head(resexp2_amce_each_municipal)
resexp2_amce_diff_municipal1 <- cj(selected ~ candgender + candparty + candage + 
                                 candexperience + candedu + candfam + 
                                 candlivetype + candorigin + candimppol, 
                               data = subset(dcj_exp2, conjoint_eleclevel=="municipal" & 
                                               genderval_2cat=="Liberal"),
                               feature_labels = vlab_exp2,
                               by = ~dv, 
                               id = ~Response.ID, estimate="amce_diff")
resexp2_amce_diff_municipal1$genderval_2cat <- "Liberal"
head(resexp2_amce_diff_municipal1)
resexp2_amce_diff_municipal2 <- cj(selected ~ candgender + candparty + candage + 
                                 candexperience + candedu + candfam + 
                                 candlivetype + candorigin + candimppol, 
                               data = subset(dcj_exp2, conjoint_eleclevel=="municipal" & 
                                               genderval_2cat=="Traditional"),
                               feature_labels = vlab_exp2,
                               by = ~dv, 
                               id = ~Response.ID, estimate="amce_diff")
resexp2_amce_diff_municipal2$genderval_2cat <- "Traditional"
head(resexp2_amce_diff_municipal2)

resexp2_amce_both_municipal <- rbind(resexp2_amce_each_municipal,
                                 resexp2_amce_diff_municipal1,
                                 resexp2_amce_diff_municipal2)
resexp2_amce_both_municipal$BY <- 
  factor(gsub("\\*\\*\\*.*$","",resexp2_amce_both_municipal$BY), 
         levels = c("Preference","Expectation",              
                    "Expectation - Preference"))
resexp2_amce_both_municipal$genderval_2cat <- 
  factor(resexp2_amce_both_municipal$genderval_2cat, 
         levels = c("Liberal","Traditional"))
resexp2_amce_both_municipal$level <- 
  factor(resexp2_amce_both_municipal$level, 
         levels = rev(levels(resexp2_amce_both_municipal$level)))

##############
## Figure 2 ##
##############

tmp1 <- subset(resexp1_both, feature=="Gender")
tmp1$exp <- "Experiment 1\n[Without Policy Attribute]\n(House of Representatives)"
tmp2a <- subset(resexp2_both_house, feature=="Gender")
tmp2a$exp <- "Experiment 2\n[With Policy Attribute]\n(House of Representatives)"
tmp2b <- subset(resexp2_both_municipal, feature=="Gender")
tmp2b$exp <- "Experiment 2\n[With Policy Attribute]\n(Municipal Council)"

fig2dt <- rbind(tmp1, tmp2a, tmp2b)
fig2dt$exp <- factor(fig2dt$exp, unique(fig2dt$exp))

setfacetfig2_vline <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             VL = rep(c(0.5,0.5,0),each=2))
setfacetfig2_lims <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             level = factor("Male", levels(fig2dt$level)), 
             feature = factor("Gender", levels(fig2dt$feature)),
             lims = c(0.45,0.55,0.45,0.55,-0.06,0.06))

## Plot!
require(ggplot2)
p_resfig2 <- 
  ggplot(fig2dt) + 
  geom_vline(data=setfacetfig2_vline, aes(xintercept=VL), 
             linetype=2, color="gray70") + 
  geom_errorbarh(aes(xmin=lower,xmax=upper,y=level, color=genderval_2cat), 
                 height=0, linewidth=0.5,
                 position=position_dodge(width=-0.8)) +
  geom_point(aes(x=estimate,y=level, shape=genderval_2cat), size=1.5,
             position=position_dodge(width=-0.8)) + 
  geom_point(data=setfacetfig2_lims, aes(x=lims,y=level), 
             size=1, alpha=0) + 
  scale_shape_discrete(name="Views on Gender Roles") + 
  scale_color_manual(name="Views on Gender Roles", values=rep("black",2)) + 
  facet_grid(exp~BY, scales = "free", 
             switch = "y") + 
  scale_x_continuous(breaks = c(-0.05,0,0.05,0.47,0.5,0.53)) + 
  labs(x="Marginal Means Estimates (with 95% Confidence Intervals)",
       y=NULL) + 
  theme_bw(base_family = "serif") + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        panel.grid = element_blank(),
        strip.placement = "outside",
        strip.text.y.left = element_text(angle=0),
        strip.background.y = element_blank())
p_resfig2
ggsave("fig2_v4.pdf", width=7, height=3.5)

#######################################
## Appendix D.2 = Figure 2 AMCE ver. ##
#######################################

tmp1 <- subset(resexp1_amce_both, feature=="Gender")
tmp1$exp <- "Experiment 1\n[Without Policy Attribute]\n(House of Representatives)"
tmp2a <- subset(resexp2_amce_both_house, feature=="Gender")
tmp2a$exp <- "Experiment 2\n[With Policy Attribute]\n(House of Representatives)"
tmp2b <- subset(resexp2_amce_both_municipal, feature=="Gender")
tmp2b$exp <- "Experiment 2\n[With Policy Attribute]\n(Municipal Council)"

fig2_amcedt <- rbind(tmp1, tmp2a, tmp2b)
fig2_amcedt$exp <- factor(fig2_amcedt$exp, unique(fig2_amcedt$exp))

setfacetfig2_amce_vline <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             VL = rep(c(0,0,0),each=2))
setfacetfig2_amce_lims <- 
  data.frame(BY = factor(rep(c("Preference","Expectation",              
                               "Expectation - Preference"), 
                             each = 2),
                         levels = 
                           c("Preference","Expectation",              
                             "Expectation - Preference")),
             level = factor("Male", levels(fig2_amcedt$level)), 
             feature = factor("Gender", levels(fig2_amcedt$feature)),
             lims = c(-0.08,0.08,-0.08,0.08,-0.14,0.02))
fig2_amcedt

## Plot!
## Figure A9 (Appendix D.2)
require(ggplot2)
p_resfig2_amce <- 
  ggplot(fig2_amcedt) + 
  geom_vline(data=setfacetfig2_amce_vline, aes(xintercept=VL), 
             linetype=2, color="gray70") + 
  geom_errorbarh(aes(xmin=lower,xmax=upper,y=level, color=genderval_2cat), 
                 height=0, linewidth=0.5,
                 position=position_dodge(width=-0.8)) +
  geom_point(aes(x=estimate,y=level, shape=genderval_2cat), size=1.5,
             position=position_dodge(width=-0.8)) + 
  geom_point(data=setfacetfig2_amce_lims, aes(x=lims,y=level), 
             size=1, alpha=0) + 
  scale_shape_discrete(name="Views on Gender Roles") + 
  scale_color_manual(name="Views on Gender Roles", values=rep("black",2)) + 
  facet_grid(exp~BY, scales = "free", 
             switch = "y") + 
  scale_x_continuous(breaks = c(-0.1,-0.05,0,0.05)) + 
  labs(x="Average Marginal Component Effect (with 95% Confidence Intervals)",
       y=NULL) + 
  theme_bw(base_family = "serif") + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        panel.grid = element_blank(),
        strip.placement = "outside",
        strip.text.y.left = element_text(angle=0),
        strip.background.y = element_blank())
p_resfig2_amce
ggsave("fig2_amce_v4.pdf", width=7, height=3.5)

### Extracting p-values
### Figure A10 (Appendix D.2)
require(ggplot2)
p_resfig2_amce_pval <- 
  ggplot(fig2_amcedt) + 
  geom_label(aes(x=ifelse(genderval_2cat=="Liberal",0.49,0.51),
                y=level,label=sprintf("%.03f",p), 
                fill=genderval_2cat), size=3) + 
  geom_point(aes(x=0.48,y=level), size=1, alpha=0) + 
  geom_point(aes(x=0.52,y=level), size=1, alpha=0) + 
  scale_fill_discrete(name="Views on Gender Roles") + 
  facet_grid(exp~BY, scales = "free", 
             switch = "y") + 
  labs(x="p-values from significance test (H0=0)",
       y=NULL) + 
  theme_minimal(base_family = "serif") + 
  theme(plot.title = element_text(hjust=0.5), 
        legend.position = "bottom",
        strip.placement = "outside",
        strip.text.y.left = element_text(angle=0),
        strip.background.y = element_blank(),
        axis.text.x = element_blank())
p_resfig2_amce_pval
ggsave("fig2_amce_pval_v4.pdf", width=7, height=3.5)
